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Abstract. We study the percolation properties of force networks in an anisotropic 
model for granular packings, the so-called q-model. Following the original recipe of 
Ostojic et al. [Nature 439 828 (2006)], we consider a percolation process in which forces 
smaller than a given threshold / arc deleted in the network. For a critical threshold 
/c, the system experiences a transition akin to percolation. We determine the point 
of this transition and its characteristic critical exponents applying a finite-size scaling 
analysis that takes explicitly into account the directed nature of the q-model. By means 
of extensive numerical simulations, wc show that this percolation transition is strongly 
affected by the anisotropic nature of the model, yielding characteristic exponents which 
are neither those found in isotropic granular systems nor those in the directed version 
of standard percolation. The differences shown by the computed exponents can be 
related to the presence of strong directed correlations and mass conservation laws in 
the model under scrutiny. 
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1. Introduction 

Granular media sliow a peculiar kinetic behavior including the possibihty of exhibiting 
jammed configurations. Jammed assembhes of grains at high densities are not able 
to explore the phase space but can eventually yield at hirfi drives, for instance under 
shear stress, like a viscoplastic solid or a complex fiuidU, Q. Over the last years, 
experimental observations and numerical simulations of jammed granular media have 
repeatedly shown the heterogeneous distribution of stress and contact forces in dense 
packings Starting from the first studies of weight distributions in bead packs 0, IHl; 
the presence of force chains has been especially emphasized, chains which form an 
intricate force network structure and are responsible for most of the material's unusual 
properties. Force networks in a dense granular packing play the role of the cytoskeleton 
in a living cell, thus determining its mechanical response and stability. They are also 
at the core of several important properties of granular media such as friction and wear 
a, =ou„d transmission g or even electrical transport g 

Internal stress and contact forces can be determined experimentally, using for 
example photo-elastic materials, which exhibit stress-induced birefringence. The results 
obtained from birefringent packings confirm that large forces seem to indeed concentrate 
along branching-like paths, i.e. force chains or arches. Following some of these 
measurements, it was argued that a close inspection of contact force properties (for 
instance, the shape of force probability distributions) could provide new insights 
regarding the jamming- yielding transition in granular matter j^, Eoj. Nevertheless, the 
distribution of forces alone does not describe the rich topological features observed in 
experiments nor their potential physical consequences, and complementary methods are 
thus required for their analysis. 

The force network in a granular system is usually defined by the contacts exerted 
between pairs of particles in the bulk of the system, in such a way that, if particles i and 
j are in contact, they mutually exert a symmetric force ftj that can also include elastic 
and/or friction interactions. We can represent these pairwise interactions in terms of 



a graph or network in which vertices represent the particles, and two vertices are 
joined by an edge if the respective particles are in contact. This force network can be 
further characterized as a weighted network, in which each edge has assigned a real value 
fij, representing the actual value of the force exerted by the vertices (the particles) i 
and j at the ends of the edge. 

Recently, Ostojic et al. 12, 13| proposed a novel way to obtain information about 
the structure of force networks in static granular matter. The method is formulated 
in analogy with percolation theory Il5| and is based on the scaling properties of 
clusters of particles connected by relatively large forces. Since each edge carries a force 
fij, a natural way to visualize the paths that carry the largest weight (arches) is to 
consider only those edges with a force larger than a given threshold /, fij > /, deleting 
those with fij < f. For small values of /, essentially all forces remain in the system, 
and they form a connected network with a single cluster encompassing all the particles 
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in the system. Upon increasing the value of /, the network is expected to break down 
in subnetworks of connected forces, each representing a path of large weight. Each one 
of these subnetworks can be understood in terms of clusters in a percolation problem 



14( 1 . By analogy with the standard percolation transition, one expects to find a critical 
percolation threshold fc, such that for f > fc the force network is fragmented into a large 
number of small clusters, while for f = fc a large spanning cluster develops, reaching the 
boundaries of the system. This analogy with a percolation transition makes it possible 
to characterize complex contact force networks in terms of a reduced number of critical 
exponents 



15 



In Ref. 121] the percolation transition in contact force networks was first studied by 



applying a finite-size scaling (FSS) [16[ data collapse technique. This technique allows 
to estimate the value of the percolation threshold fc as well as some exponents related 
to the divergence of the average cluster size in the infinite network size limit. The 
remarkable conclusion of this work is that different isotropic models of a dense granular 
packing seem to exhibit similar percolation exponent values, independently of their 
microscopic details. Thus these exponents appear to define a robust new universality 
class for contact force networks, a class which, on the other hand, is different from that 
of standard percolation. 

Many real granular systems, however, are strongly anisotropic; for example, sand 
piles and silos are driven by the action of gravity, and have therefore a preferred 
(downwards) direction. The presence of anisotropy should in these other cases be 
naturally reflected in the contact force network percolation transition, making it in 
principle more akin to the anisotropic counter par t of percolation, namely directed 



percolation 17|. In fact, in Ref. 13|] (see also 12|) it was observed that anisotropic 



packing models indeed exhibit a different scaling in their force network percolation 
transitior(||. The results in Refs. 12, 13| however, were based in the application of an 
intrinsically isotropic formalism to an anisotropic system, not taking into account, for 
example, that correlation lengths along different directions might scale differently. 

Our purpose in this paper is to fill in this gap, presenting a detailed study of 
the force network percolation transition in an anisotropic system, performing a direct 



anisotropic scaling analysis. We focus on the q- model [19[ , a toy granular model intended 
to represent the behavior of silos, having a clearly defined preferred direction, in which 
the weight of the particles is transmitted by virtue of gravity. Performing a detailed 
FSS numerical analysis we uncover the anisotropic nature of this model, which shows 
up mainly in the presence of two correlation lengths, with different scaling behavior 
near the percolation threshold. Our numerical simulations allow us to determine a 
number of critical exponents, which we compare with those of directed percolation. The 
quantitative differences observed in the exponents clearly indicate that the contact force 
network percolation transition in granular systems with a preferred direction belongs to 

I On the contrary, Ref. [l^ considered granular packings under the anisotropic effects induced by the 
application of a shear stress, concluding that shear-induced anisotropy was not enough to modify the 
universal exponents observed in the isotropic case. 
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a new anisotropic universality class, which we fully characterize in terms of its critical 
exponents. 

The present paper is organized as follows: In Sec. [2] we briefly review the definition 
of the q-model used in our study. Section [3] describes the main elements of the FSS 
theory for anisotropic systems. The results of our analysis are presented in Section HI 
Finally, in Sec.[5]we summarize our results and present our conclusions and perspectives. 



2. The q- model 



The q-model [19| is defined on a tilted two-dimensional square lattice, whose sites are 
labeled by two integer numbers, {x\\,x±), x\\ = 1, . . . L\\, Xj_ = 1, . . . L±, giving its vertical 
and horizontal position, respectively. Each site in the row x\\ supports the weight of its 
two nearest neighbors in the immediate upper row x\\ — 1. Simultaneously, its own total 
weight is distributed between its two nearest downward neighbors located in row xy -|- 1. 
The transmission of weight from one row to the next is thus given by the equation 

+1 

w{x\\,x±) =wo + + lo^i^W - 1, 2;_L - a)w{x\\ - 1, x_l - a) (1) 

a=~l 

where wo is the constant weight contributed by each single site, P is a constant 
pressure downwards applied at the topmost row, and qa{x\\,x±), with a = ±1, are 
uniformly distributed random numbers between zero and one, restricted by the mass 
conservation condition J2a(la{x\\,x±) = 1. Eq. ([Tj) determines the set of weights 
w{x\\,x±) corresponding to an equilibrium configuration, as well as the corresponding 
force network. For instance, the relative forces between a particle at {x\\,xj_) and its 
upward neighbors (xy — l,xj_ — a) are given by gQ,(a;|| — l,x± — a)w{x\\ — l,x± — a). 
In the following we will consider the q-model defined on a lattice with periodic 



boundary conditions along the x± axis [20| , with massless particles Wq = and constant 
P. In this system of linear dimensions Ly and L± contains L\\L±/2 particles, 

each row bears an average constant weight per particle P, i.e. no weight is lost at the 
system boundaries, and the average force between particles is {fij) = P/2. Obviously, 
the pressure P is just a rescaling factor in all forces, so we set it equal to one, without 
loss of generality. 



3. Anisotropic finite-size scaling analysis 

In this section we review the FSS theory needed to analyze the force network percolation 
transition in an anisotropic system such as the scalar q-model. Let is first consider the 
isotropic case, in which there is a single correlation length ^, diverging as ^ ~ A^" as a 
function of the distance to the percolation threshold A = \ f — fd- Information about 
the position of the critical point and exponent values can be obtained by studying the 
normalized cluster number n{s, /), defined as the number of clusters of size s per lattice 
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site 15|. For this purpose, we define the average cluster size (or susceptibihty) 

x{f) = T.MsJ). (2) 

s 

In an infinite system, and close to the percolation threshold, the susceptibility diverges 



^ xif) ~ A~^. In a finite system of length L, the FSS hypothesis [16| states that 
the only relevant length scale is ^, and that the system size dependence can only enter 
through the ratio ^/L. Thus, at finite L the susceptibility scales as 

x{f,L)=L^/''Xo{A^L), (3) 

where Xoi^) — )■ x~^l^ for x — oo, and Xo(3;) — >■ const, for x — 0. Thus, for A = 0, 
V) would grow as a pure power law with L, while for A 7^ it would deviate 
from the power law behavior and saturate to a constant value for sufficiently large 
L. An estimate of fc can be obtained as the the one yielding the best power law fit 
to ^-s a function of L. Once /c is determined, a linear regression provides an 

estimate of the exponent ratio 7///. Additional exponents (and exponent relations) can 
be computed from a closer examination of the normalized cluster number. In fact, close 



to the percolation threshold, the normalized cluster number scales as 

n(s,/) = s--^(sAi/'^), (4) 

where a is a critical exponent giving the characteristic cluster size, s^ ~ A~i/'^, and T 
is a universal function, independent of s and A. Substituting A ~ ■, and defining 
the fractal dimension D as Sc ~ one obtains D = l/i^au). Right at the percolation 
threshold, in a system of finite size L, the cluster number will scale as 

n{sJ,,L) = s-^f{sL-''), (5) 

and, from Eq. ([2]), and comparing with Eq. ([3]), we obtain the scaling relation 

7 = = (3 - t)Du. (6) 

a 

The critical point and some critical exponents can also be estimated by means of 



a bisection method 2l|, |22| . Consider a random realization of a force network with size 
L and an initial guess for the percolation threshold = fmaxf^, where fmax is the 
maximum force present in the network. We can estimate the true percolation threshold 
by an iterative procedure. In any step with a guess value we check whether a 
percolating (spanning) cluster exists or not. If it does, we increase the threshold by 
/c"^^ = fc + /max2"(*+^^; otherwise we decrease it as f'+'^ = fl - /max2"(*+^). Iterating 
this scheme a sufficient number of times, we compute the percolation threshold for a 
given network realization. Averaging over many random networks, we can obtain an 
estimate of the threshold {fc{L)) for the system size considered. The fluctuations of 
this estimate, cr(L) = [(/c(L)^) — {fc^L))^]^^"^, as a function of L, yields the value of the 
correlation exponent, 

a(L) ~ (7) 
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while the percolation threshold can be obtained from the average value as 

(8) 

In anisotropic systems, the FSS theory takes a slightly more complex form. The 
length of a typical cluster is now given by the correlation lengths along the longitudinal 
(downwards) and transverse directions, and respectively, that scale as 

eii-A-'^ii, and a-A-*^^, (9) 

where the exponents and are, in principle, different. The anisotropy exponent, 
measuring the differente scaling of both correlation lengths, is defined as the ratio 

(10) 



e = X 



In finite size simulations, two different length scales are thus present, Ly and Lj_. Varying 
them independentl y w ould lead to an uncontrolled scaling of the relevant functions. 
A proper analysis 23| shows, however, that when the longitudinal and perpendicular 
lengths are related by the constraint 



22, 



24 



25 



^11 



^±5 



the system behaves as if effectively isotropic, and standard FSS applies in terms of a 

single length scale. This fact suggest an efficient way to compute the critical percolation 

exponents by performing numerical simulations for systems with freely varing Ln, and 
1 /f) 

fixed = L\\ ■ With now a single characteristic length, the percolation threshold and 
the exponent ratio 7/z^|| can be found by a standard FSS analysis of the susceptibility. 



which at the critical point takes the form 23 



Analogously, the normalized cluster number will take the form 



n{sJ,,L\\,Lf)) = s-^f{sL. 



-D« 



:i2) 



(13) 



where the exponent D\\ will satisfy the anisotropic equivalent of Eq. ([6]), namely 

7 = (3 -r)D||i/||. (14) 

The bisection method described above can also be analogously modified to work 
in anisotropic systems 22|. With the rescaling of system lengths given by Eq. (fTTl) the 
variance of the threshold estimate at finite sizes takes the form 



a{L\\,L\{") = [{f,m,L\/y) - (fcm^Lf)) 



1/2 



-l/^ll 



and the percolation threshold is given by 



|A-(/,(L||,Ln)|~L, 



(15) 



(16) 
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Figure 1. a) Correlation lengths in the q-model, computed at fixed = 16384 
and variable L± ^ L||, for different values of the threshold force /. b) Susceptibility 

1/9 

xif, , L|| ) of the q-model for different values of /. The straight line corresponds to 
the best power-law fit, corresponding to / = 0.690 and yielding a slope j/i'w — 1.37. 



4. Computer simulations 

We have studied the percolation transition in the force network of the scalar q-model 
by means of computer simulations on systems of size Ly x L±, with Ly up to 66451 and 
L± up to 2048. In order to apply the anisotropic FSS scheme described above, they 
key point is to have an a priori knowledge of the anisotropy exponent 6. A numerical 
estimate of this exponent can be obtained using the fact that, close to the critical point, 
the two correlation lengths must be related by ~ ^l- Consider a system of very 
large longitudinal size Ly and a small transversal size L± <^ Ly. In the vicinity of the 
percolation threshold, for small Lj_ we will have ^± ~ L±, and by increasing Lj_, we 
will observe that ^y increases as ^y ~ ~ For sufficiently large Lj_, and not too 
close to the threshold, we will have that both and ^_|_ saturate to their corresponding 
values given by Eq. (j9]). Therefore, the exponent 6 can be determined by simulations 
at fixed and large Ly, by plotting as a function of ^± computed for increasing, but 
small, L± values, and different force thresholds. The / values yielding the best power 
law fits are in the vicinity of the percolation threshold fc- 

In Fig. [U^a), we present the results of simulations of the q-model with fixed 
Ly = 16384 and L± running from 16 up to 2048, for different values of f. Correlation 



lengths were computed as is customarily done in anisotropic systems [26|: For each 
cluster c of connected forces that is composed by a set of vertices {x^l\x^^}, with 
z = 1, . . . , s, we define the quantities 

^ii(^) = -hxT-x\;^i Rlic) = ]t{x^-x?r, (17) 

^ 1=1 ^ 1=1 

where Xy" and are the coordinates of some reference point within the cluster. We 
have chosen the point with the highest longitudinal coordinate and the average x± 
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coordinate, respectively. Then, the correlation lengths are defined as 

, j:'sR\\is)sMsJ) ,2 nRlis)sMsJ) .... 

where the prime indicates that one has to exclude the spanning clusters from the sum 
over cluster sizes. From the plots in Fig. [T](a), in which we have represented the data 
providing a best power-law fitting, we conclude that the percolation threshold is located 
in the vicinity of / ~ 0.70. Moreover, a linear regression for the smallest values of 
yields an estimate of the anisotropy exponent 6 = 1.78 ± 0.05. 

Once the exponent 6 has been estimated, we can proceed with the full FSS analysis. 
In the first place, we focus on the behavior of the susceptibility xifi computed 
for L|| ranging from 12 to 66451. In Fig. [U^b) we represent the susceptibility as a 
function of Lu for different values of /. As can be seen in the plot, the best power 
law behavior for %(/, ) is obtained for the threshold force fc = 0.690 ± 0.005; 

significant deviations can be observed for slightly larger and smaller values of /. A 

1 /6 

linear regression of ) at the percolation threshold yields the exponent ratio 

7/i/|| = 1.37 ±0.01. 

The numerical analysis of the full normalized cluster size distribution at the 
percolation threshold can be performed using the moment analysis technique developed 



for the study of self-organized critical systems 27| . The k-th moment of the cluster 
distribution is defined as 

Mfc = ^/n(s,/e,L||,L;/')), (19) 



At the percolation threshold, when the cluster number is given by Eq. (1131) . we have 
that Mfc(L||) ~ L"^^''\ where the fc-dependent exponent is given by 

a{k) = D\\k + D\\{1 -t). (20) 

Thus, computing Mfc(L||) as a function of Ly for different system sizes provides 
information on a{k), which should be a linear function of k of the form a{k) = ao + kai, 
from which we obtain = cti and r = 1—ao/ai. The correctedness of exponent's values 
can be checked by means of a data collapse technique: Noticing that the normalized 

11 



rD\\ 

cluster number n{s, fc) scales as given by Eq.(fT3l). then L,, n{s, fc) should collapse onto 



a universal function when plotted as a function of the rescaled variable sLy ^" . 

In Fig. [2t^a) we plot the a{k) evaluated from linear regressions of the moments 
Mfc(L||), computed from numerical simulations at the percolation threshold with system 
sizes L|| = 139, 478, 1641, 5634, and 19349. A linear regression of this function, provides 
the values = ai = 1.48 ± 0.01, and ao = —1.57 ± 0.01, from which we obtain 
r = 2.06 ± 0.02. This last value can be checked against the scaling relation Eq. ([6]) 
(properly redefined for anisotropic systems), which leads to r = 3 — 7/(i^j|Z^||) 2.07, in 
perfect agreement with the estimate from the regression of the a{k) function. In order 
to check the accuracy of these exponents for the q-model we perform a data collapse 
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Figure 2. a) Plot of the a{k) functions for the q-model at the percolation threshold. 
The straight line is a least-squares fitting yielding the corresponding Z)||. b) Data 
collapse analysis of the integrated cluster number for the q-model at the percolation 
threshold. Systems sizes are L|| = 478, 1641, 5634, and 19349. 




Figure 3. a) Fluctuations of the percolation threshold estimated by means of the 
bisection method, b) Extrapolation of the critical point from the bisection method. 



analysis of the integrated cluster number at the percolation threshold, defined as 

oo 

n,(s,L||) = ^n(s,L||). (21) 

s'=s 

In Fig|2t^b) we observe that, as expected, the plots of the integrated cluster number, 
under the rescaling s —j- sLy and n^s, L\\) — )■ Ly^ ^''^"?7,c(s, Ly), collapse onto a single 
universal function for different values of Ly. 

We turn now our attention to the application of the bisection method. Fig. [3]^a) 
shows the fluctuations cr(L||, Ly^^) computed as a function of Ly. A linear regression 
provides the slope l/z^y, from which we estimate the corresponding critical exponent 
z/|| = 3.77 ±0.01. With this result, we can compute the exponent 7 from the ratio 7/^^11, 
obtaining 7 = 5.18 ± 0.04, and from Eq. f fTOj) . u± = 2.12 ± 0.01. Finally, using the 
previously computed exponent, we can plot (/c(Iv||, Ly^^)) as a function of Ly ^ as in 
Fig. El^b), which shows a good linear behavior, with an intercept with the vertical axis 
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Table 1. Critical exponents for the percolation transition in the contact force network 
of q-modcl, compared with the values corresponding to isotropic percolation (IP) in 
two dimensions, directed percolation (DP) percolation in 1+1 dimensions, and the 
percolation transition in isotropic contact force networks (ICFN). Exponents from 

Refs. EliiEi- 



Exponent 


q- mo del 


IP 


DP 


ICFN 


7 


5.18 ±0.04 


43/18 = 2.3889 


0.54386 


2.8 ±0.2 


Z/|| 


3.77 ±0.01 


4/3 = 1.3333 


1.73383 


1.6 ±0.1 




2.12 ±0.01 


4/3 = 1.3333 


1.09684 


1.6 ±0.1 


D\\ 


1.48 ±0.01 


91/48 = 1.8958 


1.4727 




T 


2.06 ±0.02 


187/91 = 2.0549 


2.108 






1.37 ±0.01 


43/24 = 1.7917 


0.3137 


1.78 ±0.02 



providing the value fc — 0.689 in excellent agreement with the threshold obtained from 
the analysis of the susceptibility. 

5. Summary and discussion 

In Table [H we summarize the results we have obtained in our percolation analysis of 
the contact force network in the anisotropic q-model, compared with the exponents 
for isotropic and directed percolation, and with the exponents (or exponent ratios) 



available for the percolation transition in isotropic contact force networks [12|, Il3| . We 



note that the results obtained here for the q-model are compatible with those reported 



in Refs.|12l. Il3|. namely fc — 0.70, 7/// ~ 1.38, and i/ ^ 3.1. Our method for estimating 
exponents is however more accurate and systematic, being at the same time capable of 
providing new exponents, not considered previously. This is specially evident for the 
exponent z/ = 3.1 ± 0.1 quoted in 12|, which is not discerning between the parallel and 
perpendicular directions. 

The main conclusion extracted from the analysis of these exponents is that, at least 
in two dimensions, the percolation transition in the contact force network of anisotropic 
granular matter belongs to a universality class different from either anisotropic contact 
force networks and isotropic percolation. It is noteworthy that the change of universality 
goes thus beyond the simple presence of a preferred direction, as we can see from the 
comparison of the q-model exponents with those of directed percolation. Even though 
some exponents are similar, such as r or D\\, others are clearly different, out of the 
estimated error bars. The ultimate reason for this difference can be traced out in the 



presence of force correlations or arches 20|. The strength of these arches is enhanced 



in anisotropic models with a preferred direction for the propagation of weight, and 
explain the change in universality between different packing models. The origin of 
correlations is easy to understand in the present case: As we have defined it, the total 
force between rows is constant, imposing a global conservation law, superimposed to 
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Figure 4. Standard deviation of the force distribution as a function of the 
longitudinal coordinate a;|| article forces for two values of the longitudinal direction 
in the masslcss q-modcl of size 500 x 500, averaged over 500 system realizations. 
Pressure per particle P = 1.. Inset: Average nearest neighbor forces as a function of 
the force /. 



the local conservation of weights built in the definition of the model, Eq. ([T]). Global 
conservation prevents dissipation of stresses, and as a consequence any local build up 
of forces will propagate downwards unchecked and lead to the creation of arches in 
which strong forces are preferably connected to one another. The strong anisotropy and 
correlations in the force network of the q-model are checked in Fig. HI where we plot the 
variance of the force distribution, computed at different heights x\\, and which shows a 
marked dependence on this variable. In the case forces were uncorrelated at different x\\ 
levels, and considering that forces are exponentially distributed 12| , the variance should 
take the form = 0.25, clearly smaller than the numerically computed values. On 
the other hand, correlations between nearest neighbors are checked in the inset of Fig. HI 
where we plot the average value {f)nn of the forces connected to a given bond of force 
/ 29|]. As we can see, this average value grows almost linearly with /, while in absence 
of correlations it should be equal to the average force (fij) = P/2 = 0.5. 

We conclude therefore that force networks in granular matter define different 
universality classes, depending on the symmetries imposed on the systems, universality 
classes that bear no resemblance with the corresponding ones in standard percolation, 
and are strongly affected by the strength of correlations in the overall force network 
structure. This result calls for further research in order to clarify the situation in more 
realistic settings, where the anisotropy might not be as strong as in the simple q-model 
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